{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Fitting step responses\n",
    "\n",
    "It is often prohibitively expensive to develop first principle models of processes and therefore it is very common to estimate low order transfer functions directly from plant data. This is simple to do if we have access to step test results."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 29,
   "metadata": {},
   "outputs": [],
   "source": [
    "import control\n",
    "import numpy\n",
    "import matplotlib.pyplot as plt\n",
    "import tbcontrol\n",
    "tbcontrol.expectversion(\"0.1.10\")\n",
    "%matplotlib inline"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Let's start with a higher order process to generate our \"real data\""
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 30,
   "metadata": {},
   "outputs": [],
   "source": [
    "Greal = control.tf([1, 2], [2, 3, 4, 1])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 31,
   "metadata": {},
   "outputs": [],
   "source": [
    "ts, ys = control.step_response(Greal)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Remember that the real data will not necessarily start at zero, so we'll add in some value for the initial output. We will also add some normally distributed noise to represent measurement error."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 32,
   "metadata": {},
   "outputs": [],
   "source": [
    "yinitial = 10\n",
    "measurement_noise = numpy.random.randn(len(ys))*0.05"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 33,
   "metadata": {},
   "outputs": [],
   "source": [
    "ym = ys + yinitial + measurement_noise"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 34,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "[<matplotlib.lines.Line2D at 0x1c1f8f10b8>]"
      ]
     },
     "execution_count": 34,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXoAAAD4CAYAAADiry33AAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjAsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+17YcXAAAgAElEQVR4nO3deXxU9b3/8dcnk21YA4JKBhBERFQq0bi0KFIUobW3BGqrtrX2V7n0tnqXLtzirW3tYqHluvRWb+vGxbZqtRUj1iUFRUXrAhhscKGgsiVA2MI6WSbz/f0xMzhJZpLJMplk5v18PPKYmTNncr4cJu/5znc75pxDRETSV1aqCyAiIsmloBcRSXMKehGRNKegFxFJcwp6EZE0l53qAsQyZMgQN2rUqFQXQ0Sk11i7du0e59zQWM/1yKAfNWoUa9asSXUxRER6DTPbEu85Nd2IiKQ5Bb2ISJpT0IuIpDkFvYhImlPQi4ikOQW9iEiaU9CLiKS5HjmOXkSkK5WWV7KobANVNX4KC7zMmz6OkiJfqovVbRT0IpLWSssruXFpBf6GRgAqa/zcuLQCIKGwT4cPCQW9iPQKHQ3cRWUbjoV8hL+hkUVlG9p8fWc/JHoKBb2IJCxVtdvOBG5Vjb9d26N15kOiJ1HQi0irIuFeWePHgMjFR7uzdtuZwC0s8FIZI9QLC7xtHrczHxI9iYJeROJqXpNufoXp7qrdJhq4sb5xzJs+rsm/AcCb42He9HFtfkPpzIdET6LhlSISV6yadHNdXbstLa9k0sLnGT3/KSYtfJ7S8sq4wRq9PfKhVFnjx9H0G8eC2RPwFXgxwFfgZcHsCQAt9v/WI+sYFXXcedPH4c3xNDlm5EOiNzHnmn9Gp15xcbHTMsWSDnr7iI3R859qUYtvzlfg5ZX5U7vkeM2/QUAoWD93jo/H1lY22R5pRvKFz2ukeSnR8k1a+HzM/aOPG/lA6A3/h2a21jlXHOs5Nd2IJEk6jNiI13QR0dW123ht8Svf282C2RNa7SuI982jo+3skWapV+ZP7TX/X/Go6UakHWI1K8TTWgdibxGr6cLCt5EmkK4MwdZCuaTIxyvzp+Ir8MbsK/CYxXxtIs0+7SlPvPdAe94b3U01epEEtbeGng4jNiL/ru5quoj3DcIRamqZN31c3PPX6BzeHE/MTtdYYnXSxipPtHjvgTVb9jVpWupp397abKM3s8XAZ4Bq59yZ4W2LgH8C6oH3gf/nnKuJ8doZwK8AD3Cfc25hIoVSG730RK216fpiBGC8/buyTbun6Kq+iFht9NG8OR7yc7LYf7Qh5vMF3hzMoOZoAwOj7scrU7yho5FjNf/GEu//1GNGY4wsTfT/uivOX2fb6JcAdwK/i9q2HLjRORcws18ANwLfa3ZQD3AXMA3YDqw2s2XOuXfaVXqRHqK1mnisGlxnhvV1tUSO19EydaQvIt6xor9BxApUf0MjedlZLWruETX+Brw5Hr50wci4NezI729+7ET+/a19m4glkW9v3dGXk9CoGzMbBfwlUqNv9tws4Arn3Jeabf84cLNzbnr48Y0AzrkFbR1PNfrM1ZU1w64Ituh9suLU2qI1r8HFOgYQ8wMg0fbuRM9Re2qr8Ua7JFKm9n5zSfRYTUb8OEd2sBFPsJFsF+SXJadz5/INVNccJcsFQ8+5IFnBIB4XJMcFIdhIlnPHtme5IANzswgEGmloCJDlHFkuiDfLuG7SSUwaPQiCwdg/zkEwyI9K11NzpJYs5zBc6NYF8QAEgxiQ5UK35oKYcwzy5nDJaUP5WOEAKrbVsPK9XRw+Wo83JwvDUVvfiOEw56juN5g/T7i01fMXT2s1+q4I+ieBR5xzf2i2/QpghnNuTvjxNcD5zrkb4hxjLjAXYOTIkeds2RL3guaSpjoTNpHXd2WwtdWMEIsBHy68vNV9WgvGyDDBqhp/zKYHiP8hATR57ZH6AA2N8f++jwVJYyOX/uxZ9u09QH6gjvyGevIa68lvqMfnNe6cfTrU1n70U1fX5OfXT68nJxggN9BATjBAdmOA3GCAnMYAnx0/BOrrQz8NDdDQwNub9xBsaCA72EhOY4DsYCPZwUbyaGRovgcCAQgE8PvryQoGyA6HdyZYN+xUSr5yG5DYeyla0oZXmtn3gQDwYKynY2yL+65zzt0D3AOhGn1nyiW9U1cuPtXWDM5EjhVvslC89ljo3LT6yISdyG+u8Tc0ee7GpRXkZxt25DAn1h6hf90R+tcfpX/dUf7281Xk+o8yw3+YvvV++tT76Vfvp09DbeinvvbY/fyGOrwNdXgDdXBzA9TVsaK1At/b+r/nX4E6Tw51nhwCnmzqPdk0eHJw2dng2Qe5uZCTc+xnb7aXhrz+BLI8NGRlE/B4aMzyEDAPV37i5NB+Hg/b99Wy8v391DkjkOUhkOUhKzub6Wf5OHPkYPB4uPW5TeytC9JoWQQti8asLBotC7I8NFgWQbMm251lhfc1XNRzYPz5hosgKyv04/GEbs0+2ha+v/y9au5dtZmdh+sZOtDL3MljmP6xQp59eye/efEDqg7V48xwWPgWnBmWlUWjC90PHnvOcAYOAzOCUSOHunL2bYeD3syuJdRJe4mL/bVgOzAi6vFwoKqjx5P0Ea/poasXn2qussbf5siNqhp/k28GsQSd444rJ3Z6Wr23vpbj/AcYfPQAg48eZJD/IIP9BynwH6Kg9hAF/kMMrD187GdA3REG1B4mO4HabW12Lkdy8jma6+VITj7+nHyO5uax3zsAf04eR3Py8fTty+cnnwpeL3e9XsWuBqM2O4/anFzqsnOp8+TSb1B/7vzaJyAvD/LzQz95eU1+Sit2cePj62N+yxgZ4wP6xlY6NOc71+ScvR3jXJ4Z9TvHXFTJfQlOsGqtI9dX4IULLmjzvAJMGzeOaTMvarF9xogRzJhxbkKTzNrS1fMTOhT04dE03wMuds4djbPbamCsmY0GKoGrgC92qJSSNlrreErG4lPNRY5X0Ccn5h/8QG9OQkPu4g07tECA2/6wiiH7qzn9yH6OP7yPyqdr+HCoMTpwCHbt4q/bqrDd1fRpqIv5+xstiwP5/djv7c+B/H7s6zOADwcXcjCvHwfz+3Iwry+H8vpyKK9P6Da3D4fzvBzO68Ph3D4czckn4Gn9T/tYc0/43+Err+TOZv/uSPNX+QtHmTd9BCXj4/dfDPTmkJ+T1eERLvBRh2bzzsjWvs21Nvyz+KTBCfePdGWotjXJrC2xRnB1ViLDKx8GpgBDgF3AjwiNsskD9oZ3e8059y9mVkhoGOWnw6/9NHAHoeGVi51ztyRSKHXGpq+22qe7ukOwNbHa8VsbuodzFDYc5udn92eKtxa2boUtW2D7dvZteJ/GrdsYfGh/i/bkIEZNv4EMHj0CTjgBTjyRTfThr3uCfGh92NdnIPu9A9jnHcC+PgM5lNcHZ7HnMhZ4c6gLBGMGcqL/3nhBkszO21j7R47R2aGJ7ZXMEUiR17a3bwfa1x8VS6c7Y7ubgj59xftaG+l46qohfpHf2da7u3n4/edDq/EdqGbU/ipG7a9iZM1ORtTsZGTNLkYcrMbbUNv0F/Tty8HjC6mgH5X9jmNH/yFU9xvMzv7HsavfcVT3HcTevgUEszwxO9ba8wHVvNM1Xs04IifL6Jef3aKW3bwm3rzDN5E1YxLtUE7kd1aFFxVrrr2dkV2lPR9i8d6vbTX/QdMx/10xxFZBLz1GMicRxfqji/fHNtB/iFP2bmPsnq1MPLyDKbafwDvvMWz/ziY18iM5+WwtOJHdQwuZPP18GDUKTjopdDtyJBQUMOkXK9sM6/YMNYyEd1tB0N4JXPGOFy3e+PRIuSLB21o7dKzZqa39znhNHamaWJboezSRD4S2KjZdSYuaSY/R2iSi5tpbu4/VnmsNDSy+7xnGVG1kfPWHjNu9hXF7tnDC4X3H9qnNzmXzYB+bjh/D0nEXsXlwIZsLCtkyaBh7+wzEm5vdpD27ubb6B1prA050iYHIufjWI+uO7RPvuAZxA7KtTuvImjGxmlKi+0rihbPHLOZopngKW2myS9VSwIkOCkhk5FZPWc9eQS/dqj3BFq/TNvr10U0Pwwfk8uNTjKmHtsCaNbB6NTMrKphZXw+EAn3jcSNYNaqIDUNOYuOQEWwaMpKdA48nEHM0cGIdY611viXy+rY6HOOdi3gdyq2FSCKd1omsGRMvnNvTLh35nd29nk5bEg3nRD4QesqHmIJekqa1ae5t/RHHqy3dvOztY52R+Q21jN/yd87btp5zKt+lqOo9+teH/8gGDIDiYvi3f4OiIigqouxwH+Yve7fFH12glWaFRJoOEu1E7uqLW8daCqCtEElkREi8dvbossYL57bapZsfI/J7EnlPdJdEwzmRD4Se8iGmoJek6Oz6HbFqS9mNAU7Z+DYXbl7HpM1vcdaOf5AbDNBoWWwYehJPnD6FN32nsWPcWTy86JrQJJcoMwGXnZ1wOCX69TqRP+ZkXNz6gL+B26+c2K4QaWvFxuhadlvlirdPWyNOEv0ATZVEwznRD4Se8CGmoJek6MwsV/iotjT08H6mfLCaqe+v4cLN5fSv99NoWVScOIb7zy3h9RFnsHb46RzK63vstQYtQj4i0XBq79frtv6Yk3Vx6/aGSPMQS2SFx/Zoa1GySLl7ukQ/6CD1tfVEKOglKTq1FvuHH3J39QvU/+kxJm5/lywcO/odx5PjJ/Pi6HN49aSPcTC/X9yXtzdIuuMPtjPno6vbeZNdw4z8/nijUnrb9VZb0xNq64lQ0EtStHu0wY4d8Oij8NBD8MYbnAnUnHYm9196LUtHnM27Q0eH1hppQ0eDJNl/sJ0ZfdGbao7Remu505HG0UtSJDTppK4OnnwSFi+GsrLQcrBFRXD11XDFFTB6NND6mO2unnSSLJ1dmVOkLRpHL92u1drcpk3w29/CkiWwdy/4fDB/Pnz5yzB+fIvf1dMm1HSEareSSqrRS7coXbuNl3/9Bz774p+YvLmcYHY2WSUlMGcOpceNZ9GKTQmPWAHVhkWaU41eukRba6TEDF2/n/Jb/ocJv72Lkr3b2NlvMLdd+CVKz/kU3752CkCbww5VGxbpHNXoJSGJrJHSpIZ98CD85jdw++2waxcVJ4zh3vNm8fS4C48toesLd0T29mYZkZ5ANXrptETWSFlUtoGS0X3hjjvgV7+Cmhq47DKuOm4Kr42Y0GLUTGtDCxNdX15E2hZ7VolIM20Fb7+6o3zuqcUcHDYCfvxjmDIFVq+GsjK2fez8mEMjCwu8cYcX9oZJNSK9hYJeEhIvePMC9cx5Yymrfnsd3375QV4bcSaz5txJ6Q/vDK01Q2jCjzfH0+R1kfHurT0nIl1DTTeSkOazM7OCjXxu/fN86+UHKTy0h5dGFbFo8leoGDYWgOqoqf2JdKaqo1UkedQZKwmLjLoZU/4KN73wf5xa/SHrhp3KLy6+lldPOqvJvqm6OpBIplJnrHSJkrwDlLy4CJ59NjRr9ZFHuP7946g8UNtiX7Wxi/QcCvoM0OmLIe/dCzffHBou2b8/3HorXH895OUxLwMWrhLp7RT0aS6RddDj7WOBADNXPwU33QQHDsC//EtoRM2QIcd+vyYzifR8Cvo0l8g66LH2mfDBW5z+2eth5wcwdWpobPyECTGPEb3yY6xrmyr0RVJLQZ/mElkHPfr+8Yf28l8vLKbknRepHDCUNxbdzbcaxlD14FYKC3a3GtydvaqUiCSHgj7NJbIOemGBl117D/HVtcv4j1ceJqcxwK8+cRW/ueAK6vbk4wh1trYV3J29qpSIJIcmTKW5RCYkLTj+AE8/8B/ctHIxbww/g2nX/S+3X/RlanPyW6wDHwnuWDp1VSkRSRrV6NNcq52le/bAvHlMXrKEoyf6+N41P+GRYUVtXskpXnB35ipKIpI8CvoM0OIyecEgLF5M/be/Q9ahQ9x7/hX8+dNf5V//aSKPPrIu7tWcIuIFd1df21REuoaCPtO89x7MnQurVvH3EWdw4+e+ycahJ8HR0LrwBX1y2H+0Ie7LWwtuDbUU6ZkU9GmizUlR9fWwcCHccgv07cuC2d/hnlMuxtlH3TT+hkbysrPw5nia1MoNcITWiG8ruJN9kW0RaT8FfRpoc1jj66/D174G77zDthkz+cbZX2Z9Y+zmlwP+Bm6/cqJq5SJpREGfBuINa/zVX/5OyUN3wG23gc/Hq//zAF+rPr7VC4gUFnhVKxdJMxpemQZijYI5q2oD990xF/77v2HOHFi/nu8eGd5qyKvjVCQ9qUbfC8Rrf49sjx4l4wk2csPfHuFf//ZH9gwcAsuXw6WXAq2PZ0+k/V1EeicFfQ8Xr/19zZZ9PLa2skkNfXjNTn69bBFFOzaw7MypZN11J5+ZPP7Y8/HGuetC3CLpTU03PVy89veHX9/WZPu0ja/x9JJ/Z8y+7fzwqu8T/N3vmoQ8JDZLVkTSj4K+h4vX3NIYvjJYdmOA+SsXc+/Sn7F50DAu/+qvOPs/v8Gisg2Mnv8UkxY+T2l5JRAagbNg9gR8BV6MUE1+wewJaq4RSXNquunh4jW3AAyoPcxdpQu5aMs6fld0OT+bOoc+/fu0OtRSI2pEMk+bNXozW2xm1Wa2Pmrb583sbTMLmlnMaxSG99tsZhVmts7MdBHYDojV3AIwomYnS3//Xc7ftp55n/p3fnjZN/B48zEj7gqSIpKZEmm6WQLMaLZtPTAbeCmB13/SOTcx3kVrpXXRzS0RRZXvUfq7b3Pc0QNcc+VP+dPHph1rhqmJs3yBVpAUyVxtBr1z7iVgX7Nt7zrnVEXsJiVFPl6ZPxUDztu2nt8/+gMO5fVl1jX/zesjJ2DAK/OnUlLki7vgmFaQFMlcye6MdcBfzWytmc1tbUczm2tma8xsze7du5NcrN5p5u63eeDRH7Gz33F84YsL2Tw41NYeHeIaWSMizSW7M3aSc67KzI4HlpvZe+FvCC045+4B7gEoLi5ua6XczFNWxq1/+CEbBxfypS/8lL19C4CWIa4VJEWkuaQGvXOuKnxbbWaPA+eRWLu+RFu9GmbPxnP6eN6//ffkv1aNtRLiGlkjItGSFvRm1hfIcs4dCt+/DPhJso6XtjZtgssvhxNOgGee4fITT+TyKakulIj0JokMr3wYeBUYZ2bbzew6M5tlZtuBjwNPmVlZeN9CM3s6/NITgJfN7C3gDeAp59yzyflnpKnqapgxI3RFqGefhRNPTHWJRKQXarNG75y7Os5Tj8fYtwr4dPj+B8BZnSpdholevGxE/1yWPvYDhlRVwfPPw6mnprp4ItJLaQmEHiKyeFlljR8HXP3kPQwpf501N/0SLrgg1cUTkV5MQd9DRC9eNnXTG3zj9T/z4MQZXHFwdJP1akRE2ktB30NEZq76DlRz21O38fbxJ/OTS0JTDyLr1SjsRaQjFPQ9RGGBF3NBbvvLrXiCjXyzZD512bnHntd6NSLSUQr6HmLe9HF8pWI5529/m59cMpctgwpb7KP1akSkI7RMcQ9RcjxcvmoJq8eczZ8mXBpzH61XIyIdoRp9T+AcXH89OY0Bzi17lDuuKtJ6NSLSZVSj7wkeewyeeAIWLYIxYygJb9Z6NSLSFcy5nrd+WHFxsVuzJkOuU+L3hyZDDR0Kb7wB2frsFZH2M7O18a77oVRJgegZsPPeKuWb27fDgw8q5EUkKZQs3SwyA9bf0EiB/yBfXvkwK8eex4H+HzXZiIh0JXXGdrPoGbA3/O0R+tb7+flF12qMvIgkjWr03STSXFMZHgs/vGYn15Q/xZ/PvISNQ0/CNEZeRJJEQd8NoptrIr6z6g8EzcPtF34J0Bh5EUkeNd10g+jmGoCR+3fw2Xdf4oGzL2fngCEYofVstHiZiCSDgr4bNF+6YO4bSwlkZXH/uSUYoSuogxYvE5HkUNB3g+hmmSFH9vP5ihU8duYl7Ot/HM1nMWjxMhHpagr6bjBv+rhjSxp8de2T5DQG+P0nrqAxzmQ1LV4mIl1JQd8NSop8LJg9gbH5Qb7y5lO8eOaFfP2fP4UvTgesOmZFpCsp6LtJSZGP5f02MKDuCJ9cfCslRb4mNf0ILV4mIl1Nwyu7SyAAd9wBU6fCuecCHFukTIuXiUgyKei7y1/+ApWVcNddTTaXFPkU7CKSVGq66S533w0+H1x+eapLIiIZRkHfHTZvhrIymDNHK1SKSLdT0HeHe+8Fs1DQi4h0MwV9sjU0wP33h5pshg9PdWlEJAMp6JPtiSdg1y74+tdTXRIRyVAK+mS7+24YORJmzEh1SUQkQynok2nrVlixAq67DjyetvcXEUkCBX0yPfxw6PbLX05tOUQkoynok+jAfUuoGHk6o+95V2vNi0jKKOiT5PlHVzBw03s8Om4yDq01LyKpo6BPksr/vZ+AZfHUaRcd26a15kUkFRT0yRAMMvXN51g1uoh9fQY2eUprzYtId1PQJ8PLL+M7tJvS06e0eEprzYtId1PQJ8NDDxHI9/Ly6ZOabNZa8yKSCm0GvZktNrNqM1sfte3zZva2mQXNrLiV184wsw1mtsnM5ndVoXu0+nr405/Inj2LH1x1Hr4CLwb4CrwsmD1BSxKLSLdLZCnFJcCdwO+itq0HZgN3x3uRmXmAu4BpwHZgtZktc8690+HS9gKv3PdnJu3bxz/XjuGdsg26kIiIpFybQe+ce8nMRjXb9i6AmbX20vOATc65D8L7/hGYCaRt0JeWV1J7/4McycnnpdFF1IWHVAIKexFJmWS20fuAbVGPt4e3pa1bn3mHSzb8jZUnF1OXnQtoSKWIpF4ygz5Wdd/F3dlsrpmtMbM1u3fvTmKxkueEd9Yx9EgNZad+vMl2DakUkVRKZtBvB0ZEPR4OVMXb2Tl3j3Ou2DlXPHTo0CQWK3k+t+UN6jzZrBxzbpPtGlIpIqmUzKBfDYw1s9FmlgtcBSxL4vFSyzk+++EbvDa6iMN5fY5t1pBKEUm1RIZXPgy8Cowzs+1mdp2ZzTKz7cDHgafMrCy8b6GZPQ3gnAsANwBlwLvAo865t5P1D0m5t96ib+VWBnzxCxpSKSI9SiKjbq6O89TjMfatAj4d9fhp4OkOl643WboUsrIouuFaXumlTU8ikp40M7arPP44XHQRKORFpIdR0HeFTZtg/XqYNSvVJRERaUFB3wX+fveDAEx+t58uMCIiPU4iSyBIK0rLKxn62BNsGjycrQUngmbDikgPoxp9J925bB3FWyt4PmrsvGbDikhPoqDvpJPfepW8xkCLSVKaDSsiPYWCvpM+s72cg7l9WD389CbbNRtWRHoKBX0HlZZXMmnBc5z37mu8PLqIgOej7g7NhhWRnkSdsR1QWl7JjUsrOHn7Pzjx8D6eH3MuRmjFNl+BV2vQi0iPoqDvgEVlG/A3NPLJ91cD8MLJ5xwL+VfmT01t4UREmlHTTQdEOlqnvr+adcPGsqfvoCbbRUR6EgV9BxQWeBl09AATq/7BypPPbbJdRKSnUdB3wLzp47hk21tk4Vg5JnRtdHXAikhPpTb6Digp8nF244fU9BnA2yeMUQesiPRoCvqOcI6Ra16Gf/oU7//ys6kujYhIq9R00xEVFbBjB0yfnuqSiIi0SUHfEWVlodvLLkttOUREEqCmmwSVlleyqGwDVTV+Hn3sj5x6yjgG+tQmLyI9n2r0CYjMhK2s8ZPXUMvHPqxg6XFnaN15EekVFPQJiMyEBbhg63ryGht4/qSJWopYRHoFNd0kIHrG60Wby6nNzuWN4WdQr5mwItILqEafgOgZr5M/fJM3hp9BXU6eZsKKSK+goE/AvOnj8OZ4GHZwN2P3buPF0WdrJqyI9BpquklAZMbrO7c8B8CGj32cBbMnaCasiPQKCvoElRT5KMnZDsOG8YfbvwZmqS6SiEhC1HSTqGAQVqyASy9VyItIr6KgT9Rbb8GePTBtWqpLIiLSLgr6RK1YEbq95JLUlkNEpJ0U9IlavhzOOAMKC1NdEhGRdlHQJ6K2FlatUrONiPRKCvpEvPxyKOwvvTTVJRERaTcFfSJWrICcHLj44lSXRESk3RT0iVi+HD7+cejXL9UlERFpNwV9W/bsgfJyNduISK+loG/Lc8+Bc+qIFZFeS0HfluXLYeBAKC5OdUlERDpEQd8a50JBP3UqZGtZIBHpndoMejNbbGbVZrY+attgM1tuZhvDt4PivLbRzNaFf5Z1ZcG7xcaNsHWrmm1EpFdLpEa/BJjRbNt84Dnn3FjgufDjWPzOuYnhn892vJgpsnx56FZBLyK9WJvtEc65l8xsVLPNM4Ep4fsPAC8A3+vCcqVMaXkli8o2UFXj54EnH+Ic30j6jhmT6mKJiHRYR9voT3DO7QAI3x4fZ798M1tjZq+ZWUlrv9DM5ob3XbN79+4OFqtzSssruXFpBZU1fjyNASZuWsdfjj+D0nVVKSmPiEhXSHZn7EjnXDHwReAOM4tbNXbO3eOcK3bOFQ8dOjTJxYptUdkG/A2NAJy14x8MqD/KCyPPYlHZhpSUR0SkK3Q06HeZ2TCA8G11rJ2cc1Xh2w8INe8UdfB43aKqxn/s/kWbywli/O2ks6is8TNp4fOUllemsHQiIh3T0aBfBlwbvn8t8ETzHcxskJnlhe8PASYB73TweN2isMB77P6Fm9fx92GncMDbH4DKGj83Lq1Q2ItIr5PI8MqHgVeBcWa23cyuAxYC08xsIzAt/BgzKzaz+8IvHQ+sMbO3gJXAQudcjw76edPH4c3x0L/uCBOrNvDyqKZfQPwNjWrGEZFeJ5FRN1fHearFpZacc2uAOeH7fwMmdKp03aykyAfAq7cvJtsFWTWqZUtTdPOOiEhvoOmezZQU+Sjpt4OjuV7e9J3W4vno5h0Rkd5ASyA05xw88wwHPjGZ7Pz8Jk95czzMmz4uRQUTEekYBX1zGzbA5s0Mu7KEBbMn4CvwYoCvwMuC2ROONe+IiPQWarpp7tlnQ7czZlAyyqdgF5FeTzX65p55Bk47DUaNSnVJRES6hII+2tGj8OKL8KlPpbokIiJdRk03fLSQ2alrX+L/6up45ZRiJqW6UCIiXSTja/TRC5lN/vBN/Nl5fHNbP82AFZG0kfFBH72Q2WGFPXYAAAe4SURBVJQP1vDqyAkccB7NgBWRtJHxQR+Z6XrS/ipG79/BCyef02S7iEhvl/FBH5npevEHawF4MRz0mgErIuki44M+spDZZRtf4/3Bw9kyqFAzYEUkrWT8qJuSIh85B2u44OcV3HPebHwFXuZNH6eJUiKSNjI+6AEu31YOwSDfvOO7fPP881NdHBGRLpXxTTcAPP44FBbCueemuiQiIl1OQe/3h9a3mTkTsnQ6RCT9KNlWrAgtfTBrVqpLIiKSFAr6xx+HgQPh4otTXRIRkaTI7KAPBGDZMvjMZyA3N9WlERFJiowN+tLySq7/5q9h716+78ZobRsRSVsZGfSRhcyK31xJnSeHx0+YwI1LKxT2IpKWMjLoF5VtIFBby8x3XmT52As4muvF39CohcxEJC1lZNBX1fiZ+v5qBvsP8uczL2myXUQk3WRk0BcWeLmi4jl29RvMqtFFTbaLiKSbjAr60vJKJi18nrrKKj75/moeP+OTNGZ5ALSQmYikrYxZ6ybSAetvaOS6d14k2wV5LNxso4XMRCSdZUzQH7uSlHNcUbGCdcPGsnHISHwFXl6ZPzXVxRMRSZq0CfrIBb6ravwUxqihRzpaz6j+gPG7N3PTtG802S4ikq7Soo0++gLfDqis8bcYFx/paL3yrb9S58nmyfGTm2wXEUlXaRH00Rf4jmg+Ln7e9HGMrD3AFyqWs2z8FA54+6sDVkQyQloEfbzml+jtJUU+Fu9eSXYwwJ2fuBJfgZcFsyeoA1ZE0l5atNEXFnipjBH2TZpldu7klKV/gGuu4cW753Rj6UREUistavSRC3xHa9Es88tfQkMD3HRTN5dORCS10qJGH2l+iTvqZtcu+O1v4UtfglNOSWFJRUS6X1oEPYTCPlZ7e2l5JYHrb2BWbR1fHHoJV5dXql1eRDJK2gR9LKXllbz00zu57dVSlpz9GV7PPo6/L60AUNiLSMZIizb6eEr/7yluefJ2VvtO55ap1wEth12KiKS7hILezBabWbWZrY/aNtjMlpvZxvDtoDivvTa8z0Yzu7arCt6m3bv52QM/YL+3P9+YdSMNnpxjT2k2rIhkkkRr9EuAGc22zQeec86NBZ4LP27CzAYDPwLOB84DfhTvA6HL7N9PxXd/zPbTJjLkaA1fn/V99vRtekjNhhWRTJJQ0DvnXgL2Nds8E3ggfP8BoCTGS6cDy51z+5xz+4HltPzA6BrBIMyZQ+OwYUy49Waq8/szZ/YPqBg2tslumg0rIpmmM52xJzjndgA453aY2fEx9vEB26Iebw9va8HM5gJzAUaOHNn+0mRlweHDLJt4GfeOn8Y7J5zcsjBajlhEMlCyR91YjG0u1o7OuXuAewCKi4tj7tOmP/6Rb89/KuYBDLQcsYhkpM6MutllZsMAwrfVMfbZDoyIejwcqOrEMdsUr/1d7fIikqk6E/TLgMgommuBJ2LsUwZcZmaDwp2wl4W3JU1CyyGIiGSQRIdXPgy8Cowzs+1mdh2wEJhmZhuBaeHHmFmxmd0H4JzbB/wUWB3++Ul4W9KUFPlYMHsCvgIvBlqlUkQynjnXsebwZCouLnZr1qxJdTFERHoNM1vrnCuO9Vxaz4wVEREFvYhI2lPQi4ikOQW9iEiaU9CLiKQ5Bb2ISJrrkcMrzWw3sKWDLx8C7OnC4vRWOg8f0bkI0Xn4SDqei5Occ0NjPdEjg74zzGxNvLGkmUTn4SM6FyE6Dx/JtHOhphsRkTSnoBcRSXPpGPT3pLoAPYTOw0d0LkJ0Hj6SUeci7droRUSkqXSs0YuISBQFvYhImkuboDezGWa2wcw2mdn8VJcnlcxss5lVmNk6M8uo9Z7NbLGZVZvZ+qhtg81suZltDN8OSmUZu0Oc83CzmVWG3xfrzOzTqSxjdzCzEWa20szeNbO3zezfw9sz6j2RFkFvZh7gLuBTwOnA1WZ2empLlXKfdM5NzKSxwmFLgBnNts0HnnPOjQWeCz9Od0toeR4Abg+/LyY6557u5jKlQgD4jnNuPHABcH04GzLqPZEWQQ+cB2xyzn3gnKsH/gjMTHGZJAWccy8Bza9iNhN4IHz/AaCkWwuVAnHOQ8Zxzu1wzr0Zvn8IeBfwkWHviXQJeh+wLerx9vC2TOWAv5rZWjObm+rC9AAnOOd2QOgPHzg+xeVJpRvM7O/hpp20bq5ozsxGAUXA62TYeyJdgt5ibMvkcaOTnHNnE2rKut7MJqe6QNIj/AYYA0wEdgC3prY43cfM+gGPAf/hnDuY6vJ0t3QJ+u3AiKjHw4GqFJUl5ZxzVeHbauBxQk1bmWyXmQ0DCN9Wp7g8KeGc2+Wca3TOBYF7yZD3hZnlEAr5B51zS8ObM+o9kS5BvxoYa2ajzSwXuApYluIypYSZ9TWz/pH7wGXA+tZflfaWAdeG718LPJHCsqRMJNjCZpEB7wszM+B+4F3n3G1RT2XUeyJtZsaGh4rdAXiAxc65W1JcpJQws5MJ1eIBsoGHMulcmNnDwBRCy9DuAn4ElAKPAiOBrcDnnXNp3VEZ5zxMIdRs44DNwNcj7dTpyswuBFYBFUAwvPm/CLXTZ8x7Im2CXkREYkuXphsREYlDQS8ikuYU9CIiaU5BLyKS5hT0IiJpTkEvIpLmFPQiImnu/wO8srXaU2W+kQAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "plt.scatter(ts, ym)\n",
    "plt.plot(ts, ys + yinitial, color='red')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 35,
   "metadata": {},
   "outputs": [],
   "source": [
    "import scipy.optimize"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We'll fit a first order plus dead time and second order plus dead time model. The `tbcontrol.responses` library contains the analytic formulae for these step responses."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 36,
   "metadata": {},
   "outputs": [],
   "source": [
    "from tbcontrol.responses import fopdt, sopdt"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "It is very important to have a good idea of the initial parameter values. Interaction makes it easy to figure out."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 37,
   "metadata": {},
   "outputs": [],
   "source": [
    "from ipywidgets import interact"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 38,
   "metadata": {},
   "outputs": [],
   "source": [
    "def resultplot(K, tau, theta, y0):\n",
    "    plt.scatter(ts, ym)\n",
    "    plt.plot(ts, fopdt(ts, K, tau, theta, y0), color='red')\n",
    "    plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 39,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "application/vnd.jupyter.widget-view+json": {
       "model_id": "bffa6e2450304400ac3be78eeabb4a4e",
       "version_major": 2,
       "version_minor": 0
      },
      "text/plain": [
       "interactive(children=(FloatSlider(value=5.5, description='K', max=10.0, min=1.0), FloatSlider(value=5.0, descr…"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "interact(resultplot, \n",
    "         K=(1., 10.), \n",
    "         tau=(0., 10.), \n",
    "         theta=(0., 10.), \n",
    "         y0=(0., 20.));"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We can use the `scipy.optimize.curve_fit` tool to do this fit just like when we did regression without time."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 26,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "[1.954094509563327, 2.8181973943514884, 0.6113717655974688, 10.031145239708668]"
      ]
     },
     "execution_count": 26,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "[K, tau, theta, y0], _ = scipy.optimize.curve_fit(fopdt, ts, ym, [2, 4, 1, 10])\n",
    "[K, tau, theta, y0]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The parameters for the second order model should be similar, with a smaller time constant and overdamped nature."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 27,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "[1.9798056401724866,\n",
       " 0.8305183954384022,\n",
       " 1.8452576915795178,\n",
       " 0.3199080792764559,\n",
       " 10.005597510469169]"
      ]
     },
     "execution_count": 27,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "[K_2, tau_2, zeta_2, theta_2, y0_2], _ = scipy.optimize.curve_fit(sopdt, ts, ym, [2, 2, 1.5, 1, 10])\n",
    "[K_2, tau_2, zeta_2, theta_2, y0_2]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 28,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "<matplotlib.legend.Legend at 0x1c1fd1c710>"
      ]
     },
     "execution_count": 28,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAlkAAAEvCAYAAAB2a9QGAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjAsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+17YcXAAAgAElEQVR4nOzdd3hUVf7H8feZkmRSILQgCQJRpChggCB2xQKKZRELiiJY17WsLCsKq+7CugqK/tBVUFlFsICNIquwiIKrLliCoIA0aUICBAIhbVJm5v7+mBCTkMCElJkkn9fzzJPMnTtzvzOZJJ8559xzjGVZiIiIiEjNsgW7ABEREZGGSCFLREREpBYoZImIiIjUAoUsERERkVqgkCUiIiJSCxSyRERERGqBI9gFVKRly5ZWhw4dgl2GiIiIyDGtXLlyv2VZrcpvD8mQ1aFDB1JSUoJdhoiIiMgxGWN2VLRd3YUiIiIitUAhS0RERKQWKGSJiIiI1IKQHJNVkaKiInbt2kV+fn6wS6nXIiIiaNu2LU6nM9iliIiINGj1JmTt2rWLmJgYOnTogDEm2OXUS5ZlkZGRwa5du0hMTAx2OSIiIg1avekuzM/Pp0WLFgpY1WCMoUWLFmoNFBERqQP1JmQBClg1QK+hiIhI3ahXISvY7HY7SUlJJZft27cD8PXXX3PGGWfQpUsXunTpwrRp00ruM27cOBISEkhKSqJbt24sWLDgiO2nnHIKgwcP5ueffwbgmmuuISkpiY4dO9K0adOS4y1fvrxMPRs2bCApKYmePXuyZcsWzj77bAC2b9/OrFmz6uAVERERkcrUmzFZocDlcrF69eoy2/bs2cPQoUOZP38+vXr1Yv/+/QwYMICEhASuuOIKAP70pz/x0EMPsX79es477zzS09PLbAd47733uOiii1izZg3z5s0D4IsvvuDZZ5/l448/rrCe+fPn87vf/Y7x48cDlISwwyFr6NChNf8iiIiISEDUklVNU6ZMYcSIEfTq1QuAli1b8swzzzBx4sQj9u3atSsOh4P9+/cfcduQIUPo379/wC1QCxcu5Pnnn+e1116jX79+AERHRwMwZswYvvrqK5KSkpg8efLxPjURkXph/qpUzpm4lMQxn3DOxKXMX5Ua7JJEALVkVYnb7SYpKQmAxMRE5s2bx7p16xg+fHiZ/ZKTk1m3bt0R9//222+x2Wy0anXE8kYA9OrViw0bNgRUy8CBA7nnnnuIjo4uaQ07bOLEiUdtARMRaSjmr0pl7Nw1uIu8AKRmuhk7dw0Ag3omBLM0kXoaskaOhHLddtWWlATPP3/UXSrqLrQsq8LB5KW3TZ48mbfffpuYmBjee++9SgefW5Z1HIWLiDRekxZvLAlYh7mLvExavFEhS4KufoasEHLaaaeRkpLC1VdfXbJt5cqVnHrqqSXXS4+9OppVq1aRnJxcK3WKiDREaZnuKm2XxmH+qlQmLd5IWqab+FgXowd0Dkrorp8h6xgtTnXpvvvuo2/fvgwePJikpCQyMjJ45JFH+Otf/1qlx5kzZw6ffvopzz33XLVriomJITs7u9qPIyJlhcofbvlNfKyL1AoCVXysKwjVVF+g7zG9FysXSl3IGvheTW3atOHtt9/mrrvuokuXLpx99tncfvvtXHXVVce87+TJk0umcHj77bdZunRppeO1qqJHjx44HA5OP/10DXwXqSGH/3CnZrqx+O0PtwZZB9foAZ1xOe1ltrmcdkYP6BykiioWyOD8QN9jei8e3dG6kOuaCcVxQMnJyVZKSkqZbevXr6dr165Bqqhh0WspUnXnTFxaYYtJQqyL/425KAgVyWGh3qpTvmUF/EFwwuDuZeoM9D2m9+LRJY75hIqSjQG2TbyiVo5pjFlpWdYR433qZ3ehiEgd09ifmlEbgWhQz4SQClXlBTo4P9D3WG29F0M9rAYqlLqQFbJERAJQ2R/upi4n50xcqjE0AQilsTJ1KdBQFGg4qI0QEYo/m+P9nRk9oHOFLYfB6ELWmCwRkQBUNPbHaTPkFno0hiZAoTRWpi5VFn7Kbw90fFltjEML9s+m/Ji1x+avOe7fmUE9E5gwuDsJsS4M/m7U8l2zdeWYLVnGmOnAlUC6ZVndirdNAq4CCoEtwG2WZWVWcN/LgBcAO/CaZVlHToMuIlIPHP4DXfqTdV6hh4N5RWX2q6gbKBTncgpGy1pj7XINtGWlovdYRT+XQPc7mvI//4paxqDin01NnwFZUSvaO9/8esS4qqr8zoRKF3Ig3YUzgJeAN0ttWwKMtSzLY4x5GhgLPFL6TsYYOzAFuBTYBXxvjFlgWdbPNVG4iEhdK/+HO3HMJxXuV1djaI5XsLqGQmmsTF2qSiiqKBxUFlaO92dV0c/fQIWDxct3h/fr0oo5K1OP+d6pynusog8hlZ2SV98C+TFDlmVZXxpjOpTb9mmpq98A11Vw1zOAXyzL2gpgjHkX+B2gkCUiIaM6LTp1NYampludgtWyFkpjZera8Yai2gjElYWa8kHrcHd4pruo5NiBtjBV5T1WleBU3wJ5TQx8vx14r4LtCcDOUtd3AX1r4HhB8+STTzJr1izsdjs2m41XX32Vvn37UlhYyMMPP8y///1vbDYbp556KlOmTKFt27YA2O12unfvjsfjoWvXrsycOZPIyMiS7UVFRTgcDoYPH87IkSNZsmQJjzzibxj85ZdfSEhIwOVy0aNHD958880yNY0ePZqFCxcycOBATj75ZCIjI7n11luZMWMG/fv3Jz4+vs5fJ5H6orr/wAINDdUJF7XxT7ayf2qpme4aHcRf0X4TBnev8L6N8cSAyp6zZYHXCz4fPP3xZnJzAcsOGCwLctyGJ+ds5cw2Cfh8VHixrMqvb9vgxMIJlim+gGX5l3trFRXO/uxCWkSHk1/oIzvfU7IPFO+L+S2NFd+2GehyyxoO5BbSLDKMAzmxQKz/dvzHB8MmC95803/98MX+y0kczC0qlfDKPX4xp83GqafH889/lr3/4cvh45T+vnlzuPPOWvnxBSSgebKKW7I+Pjwmq9T2R4FkYLBV7oGMMdcDAyzLurP4+jDgDMuyHqjkGHcDdwO0a9eu944dO8rcHuy5nVasWMGoUaP44osvCA8PZ//+/RQWFhIfH89DDz3EwYMHmTZtGna7nTfeeIOXX36Zb7/9FmMM0dHR5OTkAHDzzTfTu3dvRo0aVWZ7eno6Q4cO5ZxzzmH8+PElx73wwgt59tlnK11up0mTJuzbt4/w8PAy2492v2C/liKhoibmG6rtGbpro8aKxpLBkS0ZFc3lFOicT4HuV3rfvEIveG1YXhsRNgeXnxbP1xsPsjezgFZRLm7tm8h5J7emqIhKLx7Pkd97PP7LjzsO8dXGDLLyvESHOUlu1wKf1/D9toPk5vuIdDg4rU0sPq9h7a4s8gp8RNgddGwVQ6soF2kH89mankdBoY8wu534ppHERoTj9fof3+v1XzJzi8jILsLjtbBjIybCSZjNUXK7zwf5hT7yC33+UOD7Lexg6Xy0mtSlC6xfX/vHqfF5sowxw/EPiL+4fMAqtgs4sdT1tkBaZY9nWdY0YBr4JyM93rpqy+7du2nZsmVJmGnZsiUAeXl5vPHGG2zbtg273X+2x2233cb06dNZunQpF198cZnHOe+88/jpp5+OePy4uDimTZtGnz59GDduXKWLSJd29dVXk5ubS9++fRk7dizr168nOjqaDh06kJKSws0334zL5WLFihW4XPWriVWkNlRnsG9lAu0GqspYm0BqCbTGilrCnDaD024o8vr/1Fo+A0U2vB47lseOVWTH57GR77Hz6JR0HIMSyM+H/HwYNy+LA1lt/ft5beCxcdBr465PPdxjpZLjtoiwOfEUOigo7OUPTR4bls//deirdk6IhsJCi8JCKCyE7NzW+LxtjggYm0t9nwqsDugZH01ToCkGH7nGwyJ82I0XO02xGS9FePjJ5ODASzgeXMaHgzwObDpIrt0HPg/t8GLDix0vjjQvEQ6IsFs4jAc7XnxeD0VFhSX72PHicHuJdBoi7D7seLHhw11UgMGL3Xix2b0l2x3Gomk4GMvCjpf8ggKwvDiMD4MPu+Xfz258NAkzGOu3x7Thw2Aduc3yYYxVcl+vz4fX68FgYbN8xbX6CLMbHDYLGxY2y0uhx4sp9bj+rz5sWCXH8W/jt21Y2PCW2t/CWL/d5rSB3fx2/fB+luXD8vkA//HtxsJuKNmH4sewKG6SK97PVn4/8D/f4u9taVHAnmq/c47XcYWs4rMGHwEusCwrr5LdvgdOMcYk4v/9uBEYelxVhoD+/fvz97//nU6dOnHJJZcwZMgQLrjgAn755RfatWtHkyZNyuyfnJzMunXryoQsj8fDokWLuOyyyyo8xkknnYTP5yM9PZ3WrVsfs6YFCxYQHR3N6tX+Pz3jxo0D4LrrruOll146aguYSE0L9e6eqgz2re64j0Bei0C7AQ+HQcsHVpEDX6Edq9BBi7Aoli2D7GzIyYHcbB85BwrJySgg92AhuYc85Gb7WLHhIAX5J0KRA6/HgdfrJNcbhscXRqEvnEJfGB7LWelz2Qtc9XrpLWVbwQ0+XLixU4CDfFpSSDg5hFNAGIWEU0AEBYRRQHjx9fB9/tvKX5wUEU4BTorKbDt8/fD3R7s48FR43YGn5GLz93n5BfqRvlzDX8ndCo7c9beHLO4qM0Ah2I2B4g/QHp9Vss/hHrHiiEDY4bBpDF4sPD6ruLeuJEbgtNtw+H57vMP7A2AzZa+X/woU+izcRT58loXN2IgIsxPusPn3Kd6vwOMjp8CDr/TzMYYIp40Cr4XPZ2GzGTzlXkOr1HOy22x4i/drEuEkMrxU7DiidtuRt5XaJ6fQS0ZuIT7LfvilxWCIjnCQV+jF4wOH3dAsMozoiOLj1MBSddURyBQOs4ELgZbGmF3A3/CfTRgOLClucfnGsqx7jDHx+KdqGFh85uH9wGL8UzhMtyxrXY1UPXIkrK7+Z5oykpKOuvB0dHQ0K1eu5KuvvmLZsmUMGTKEiRMn0rNnzwpbnSzLKtnudrtJSkoC/C1Zd9xxR6XHCcVljkSOJRQnMiwv0MG+1R2IXdFr8bdZKylYv5teTgeZqblkpuUx7/Nf6XHIh1XgwFMQhqcojAKPi3/+3wZetKeR5XGR7YnkkPdUsnzR5BFd5jhpwEVTS2+xARFABGEUEEVu8cVBM3KJJI+o4q+HLy7cZb4/fIkkj3DyceEmgnwiir+GmwLCTAHhppAwU4DTFGEzXnw2Gz5jw2ts+Gz+r16bvdR1Ox6bf5vN6eS0ds3BbgenE+x2lv96iByPwWezUWS347U5KLLZ8NoceGw2imwOCmx2PMaG1+7ivkt7gsPhv7/DAWFhv31f/Ljf7crmo3Xp7Hd7aBITxd48Dx67Ha/NhsfY8dgc/mMZu79Wm40iuwOPzY7v8D7Gjsfhr8GH/7lYxmCZqnfplV/S5YIAu4HtwL+r8AEm0A87YcWXowkHFlXyeDGl9qurZX4GVHKcQLq5gyWQswtvqmDz6xVsw7KsNGBgqesLgYXHXV2IsdvtXHjhhVx44YV0796dmTNncv3117Njxw6ys7OJifntbffDDz+ULBLtcrlKWpuOZuvWrdjtduLi4mrtOYjUhsrOJBq3YF3ItG5V1r1m4f9ncESNlgWZmfDrr5CaCrt3U5Sazv5d+WTsLiQjAzIy7WTkhHEgz8WBgigOFMWQXhhDnBXLQZqRSTMO0owdxHDjUWoz+GhCFk3IoimHaMohWrKfk8imickmypZDFDm4TA5R9jxinIXYrBzCbfk4HQU4nEWEOYrolBBO4gnhOCOdEBkJERHMXrufdK+DfGcYbmc4ec4I8pwRuJo1ZdKtZ0JMDERHs3h7No8v/ZWDOCiyO8GYSlv6Ktt+LIf/+Z1W7j2QXsH4rcqOkRDr4r5j/OOevyqVsZvX4D6xQ7VrPsxuDL5qfAiuaOLRQE+GCLRLurIPOyk7DrBsw77j+j0M5Nh1ddbo0X6HSwv2HHSl1c9ldY7S4lRbNm7ciM1m45RTTgFg9erVtG/fnqioKIYPH86oUaN45ZVXsNvtvPnmm+Tl5XHRRYEn+H379nHPPfdw//33BzQe61hiYmLIzs6u9uNI43K8XX6V/fHLdBeVOf07aK1b+/dzadZWInbuID5rPyfkZNAyN5OYvDxsBRFEO5qRnhtFekFTtj4WyxhvC/bRknTi2E9L9tGJ/ZzNIWIrPUQYBbQwB2jKQWJtmbS27aajfRORjlzCHXmEh7k5q1M0sc0Msc1tLNyaxkG7j/wYOzkxYWRGxpAZEUNU6zj++/iAoz6dqrQcuFal8kolA9Ap9XMYcCq427QNaMxa+WBa2UD6WJeTqHBHQHNDQdl5pMrPx3S47n5dWh3zDMiqTFGAoWRsWmXbXE471/ZOOKKeQFVn4tGqqOzDTulpF2rj97A2nktFjvaeLC9U5tOqnyErCHJycnjggQfIzMzE4XDQsWNHpk2bBsCECRN46KGH6NSpEzabjS5dujBv3rxjhqXD3YiHp3AYNmwYo0aNqpF6R4wYwT333KOB7xKw6nT5BfrHr8ZatwoLYcMG/2XLFti2zd/alJ6OlXGA7EM+0vJiSStsyW5fHLtpQxNOYA9JrKMNeziBvbQmg5YVPryTQprbDxBjz6RJWDanNt9HhxMyaNXa0KJNOC3bR9HipKa06NiMFq0dNG8OkZHhGNOGcyauL/VaRBZf/KFkaqkAlFHJ2XePXVnmJO4KVWUwfHUmwgw0zFV2JuG4q0+rVotJcvvmRw1elb1Hq9JqCUe+NhVtG9Qz4Yh6aiJc1mQQCWZLT13MsF5Ri1ltjausKQFN4VDXkpOTrZSUlDLbNO1AzdFrKRWpzriKiv7JBqqkVQWYtGg9th3bOS9nJ7c0zePUnL2wYwekp8OBA/iyckgvaMpOK4FdtC25pOK/nkY8acSTW278EoDLXkDLyCxszmyIcBPWwkbyGXGc17sFrVtDXBy0bg3f7U7jiU9/It9z7KkHAn0tjjZ1QbCmdQhEXTyXqgj0eYfi61MXKnveFSk/Rqy+KP8+q6zFs65/BjU+hYOINCxVaR0JZJLJij7lOz2F9Ni9mdP3bKbTvh20O7SH1tkHaPFcFpGFbi73GlJpy3Y68A0deI+T2c7F/Eo7dpp27LTaUkjZOeGcdi9tWxWQcKKdXolhXJlgiI+n5NKmjf8SExOOMa2Ao59tNGLOhjIBC6q+Zhoc/xIqgairMTB18VyqItD3aCi+PnWhPrb0VFUgLZ6hdGazQpZIIxBIK0OgS79U1q04YXB3fyvBr7/CV1+xedF/2bPiB+IP7qFlXiZRhW7slo9DNOUXOvILHfmCZLZwMptNR7ZyEruJx4e95FjG+EiIh/YdbPQ5Ea5tByee+NulbVto2dKOzRZZY69VTawzWNuBo7J/7kBAM7ZX9Vih8g8r0PdoXYafUHp9qjK2LVhLGdVGi2co/QzKU8gSaeACHWsV6Kf/SYs3YuXlcuGONZyxax1d922j/cHdxE3KhEI3WFbxDEkdSacL/+UyfjZd2WTrxBbrZDJ9zcs8nj0qH0dsHo7YPGKabsXRNA9HUzf2pnk4Y/LZPmkgdam+LGJc/h9LfZhGo7pq44y8hiaUW3oaw3u0PIUskQYu0IVaK/oU/Mj57bh650qY8QysWgVbt7J0z17CvB4MkE84G+jC/+jHatvp/BLVi599XdjqboOv1OzdbRPglFPgjGa5rDq0AZrk4GiWi6Opm6goiHDaKhxAnBCEBZTr6yLGwVr0uS6FWvdcfREqgbMxvEfLU8gSaeAC7v7avJlBS+Yx6Kv/+hf72r3bv44K/jEdeziB1fZkltuT+dbZh3VWN3YXtMUqnqXZ4KNrOxtJp8JNXf1rhnXpAp06QXTJOPQo5q+KYdLitCO6uUJlAeX6+o+8Jro564NQCQxSdY3lPVqaQpZIA1dR99fJ+3/lpu3fQL9n/YFq/37/yrXFdpt4UqIHk9K8HysdZ5By6BT2HnKBF/CCo6kbZ6tDNGm1BWerbGJOyOPpOxK5/ozjW8PvsONZaDmv0FPjn46D+Y/8eFvl6ks3pzRejfE9qpBVBbt27eK+++7j559/xufzceWVVzJp0iTCwsouTpCWlsYf//hHPvzww6M+3sCBA5k1axaxsZVPcFiZcePGER0dzUMPPVTl+0rj8pe+rVn59FTO3fwdXdO30Sr3IA7rt9XI8p0x/NDyalY0G8gKx3l8s+8kUvc6IRtsudC1K1zWD3r1gtNPhx494L/bDwRl8GpFrVaVqY+fjqvTKldfuzml8WiM71GFrABZlsXgwYP5wx/+wEcffYTX6+Xuu+/m0UcfZdKkSSX7eTwe4uPjjxmwABYubDArDkko+eknmDkTli2DTZu4IjeXw7PheIydrTGJfNP2Rta2uYYvD3bnhzVhFO0F9kJiIpx/EfTtC8nJ/iU9o6KOPMSgZsFp6aloTEdl6uOn4+qMWamv3ZzSeDTG96hCVoCWLl1KREQEt912G+Bfx3Dy5MkkJiaSmJjIsmXLyM/PJzc3l+nTp3PllVeydu1a8vLyGDFiBBs2bKBr165s376dKVOmkJycTIcOHUhJSSEnJ4fLL7+cc889l+XLl5OQkMBHH32Ey+XiX//6F9OmTaOwsJCOHTvy1ltvERlZc6erS/1w1C6k776DN97wh6qtW6Go1ADyJk3ISr6I/yaO4PPwgfx3bQt+/BGsnyHsF+jTB/70JzjrLDjzTDjhhOA8v0AF2jpVXz8dV3fMisYrSahrbO9RhawArVu3jt69e5fZ1qRJE9q1a4fH42HFihX89NNPNG/enO3bt5fsM3XqVJo1a8ZPP/3E2rVrSUpKqvDxN2/ezOzZs/nXv/7FDTfcwJw5c7jlllsYPHgwd911FwCPPfYYr7/+Og888ECtPU8JPeW7kGzbtpBx71Sy0tfQZOe2sqEqLo6i3meyovvdfGYu5bMvw/juO/CmgMsFZ58N48bB+ef7W6sCWW2pLmbyDvTYsZHOai1jEuoa45gVkYasXoaskSNh9eqafcykpKOvO21ZVoVrER7efumll9K8efMjbv/666958MEHAejWrRs9evSo8PETExNLAljv3r1LgtratWt57LHHyMzMJCcnhwEDjr5wrDQ8kz9ey6DvPuaadcvolr4FV1HBb7M4x8XBWWexZ8BwFjmuYuGnDj79FLIWgc3mb6kaMwYuucTfWhVearL0QMJTMOe1qejYTpvBaTdHLN5blTXyQlljHLMi0pDVy5AVDKeddhpz5swpsy0rK4udO3dit9uJqmjgCv4QFojwUv/97HY7brf/0+yIESOYP38+p59+OjNmzOCLL744vicgIanSoLNpE0yeDIsWsWzHDmz4Q1VOmItvTuzOvzufy+w2gxlx8iX8+9+w8l7/4yUkwJAhcPnl0K8fVHZORaDhKZjz2lR07CKfFXCrVTBb4I5XYxyzItKQ1cuQdbQWp9py8cUXM2bMGN58801uvfVWvF4vf/7znxkxYsRRx0ide+65vP/++/Tr14+ff/6ZNWvWVOm42dnZtGnThqKiIt555x0SEvTHtqEoH3ROWJNC1IzHKExdR1hOln8nY0hrHs9/TkrmjZ6/Y0tuV/I2tcb9/Ql4DkUy3vi7ACdMgIEDoXt3qKDB9QiBhqdgzmtT2TEOuYtY/bf+R71vfZ5ZurGNWRFpyOplyAoGYwzz5s3j3nvv5YknnsDn8zFw4ECeeuopZs+eXen97r33XoYPH06PHj3o2bMnPXr0oGnTpgEf94knnqBv3760b9+e7t27k52dXRNPR0LApMUbSd74Hfd8O4deaRtweQoBKLA74dxz4c478d54MzPeOshzr+SR9c4J+PLCwe4lKjGDAbccYG/sNlI9WfwbF128ndm6OrBWkEDDUzDHCFXn2I1xZmkRCT0m0O6supScnGylpKSU2bZ+/Xq6du0apIqOn9frpaioiIiICLZs2cLFF1/Mpk2bjphbqy7V19eywVi2DCZOxL30i5JglRPm4tu23Xg9+WqWd+jJ7EFXMHs2fPihf+L1sAgfMR33YztpFyf1yuaS05sfseir02bAcMR4pQmDux8RLM6ZuLTCAJMQ6/Iv8lysfIvQ0R6zplXn2IljPqGiv2wG2DbxigpuERE5fsaYlZZlJZffrpasWpaXl0e/fv0oKirCsixefvnloAYsCZJNm+Cxx2DhQsjNBcAKc7H4lDN58awhrG1zCkUHIsn9OYGCz9py9jP+QeoDB8INN8CVV9qIjo4D4gB/SKpovFJ5lbXeBDrAOphjhKpzbJ2lJyKhQCGrlsXExFC+VU4aiexseOop/8Sgu3f7t0VE+JPT2LF8GpXIw7PWk/FTa3I+a0thWjPAonufAkY9B4MHQ5MmFT90VcZEVbRvVQJMMMcIHe+xdZaeiIQChSyRmrZ4sb/VauVKsCz/XAq9esEjj8D11+OzDMuWwQcvw4458RQWGJwts2h/2S889mAUd17W5piHqKylprJ9K9KQB1jrLD0RCQUKWSI1ISvLH6zeegsyM/3b4uPhrrv84crlIjUVXn8Cpk+HHTv80yvcdafh9tuhZ88mGFNJs1UFKmqpqWxMVmNtvWnIIVJE6geFLJHqWLkSHnwQli/3t1o5HHDppfDMM5CUhNcLn34Kr74KH38MXq9/YtCJE2HQIH/v4fGorKWmom0KGiIiwaGQJVJVlgUzZsD48f4mKfAv+vfAAzB6NDidpKfDa0/BtGn+XeLi/DfddRecdNLRHz7QSTQra6lRqBIRCQ0KWVVgt9vp3r07RUVFOBwOhg8fzsiRI7HZbJXeZ/v27SxfvpyhQ4fWYaVSKzwef5fglCmQk+Pf1rMnPPssXHQRlgXffOO/+b33LTxFhoj2++k0dDf/eLA5TicMez90l7EREZGapZBVBS6Xi9XFiyamp6czdOhQDh06xPjx4yu9z/bt25k1a5ZCVn2WmwsjR+KdORN7URFFNjtfn3Ye+c9N5vIBvXn/m1T+ct16dn4dT+HepoS5vMQk7cJ1+jacLXIpAP4yf2eZ8VKhuIyNiIjUrAYbsmp73bK4uIGCiBAAACAASURBVDimTZtGnz59GDduHDt27GDYsGHkFs+B9NJLL3H22WczZswY1q9fT1JSEsOHD+eaa66pcD8JQRkZcPfd8NFH4PVS5AjjzV5XMuHC2yh0huNYlMUJr6XzzcJWePMScLbIpnn/NUSfmooJP745rIK5jI2IiNSsBhmy6qrL5aSTTsLn85Genk5cXBxLliwhIiKCzZs3c9NNN5GSksLEiRN59tln+fjjjwH/5KQV7SchJCsL7rgD5s4Fnw9iYnjtjGt4ste1WDY7BWmxZKV0IG9jG7b4DK6Oe4lJ3k5Eu4yA1g0sLZSWsRERkZrVIENWXXa5HF6WqKioiPvvv5/Vq1djt9vZtGlThfsHup/UncOtngfTD/J/y16h/09Lsfl85IRH8ey5Q1ly8RB2Hcwnb/MJZH+fSEFqc0xYETG9thPTawfOZnnHfezy4UmTaIqINBwNMmTVVZfL1q1bsdvtxMXFMX78eFq3bs2PP/6Iz+cjopJz8ydPnhzQflI35q9K5dEPVjFm0VRu/PFTnJaXXGcEUy64ial9BuMrdJDz+QlkpyTiORSJIzaXZhevI7r7TmzhXuzG4K1gkTwDZdbOC3QOK02iKSLScDTIkFUXXS779u3jnnvu4f7778cYw6FDh2jbti02m42ZM2fi9fpbImJiYsjOzi65X2X7SXBsf/xJvl88nUhPAW5HGC/2HcI/z74RT04k2V90IHt1O6xCJ+FtD9Dsop9xddyLKT6Z1OW0c23vhCMWaj68fdmGfcc1h5Um0RQRaRgaZMiqrS4Xt9tNUlJSyRQOw4YNY9SoUQDce++9XHvttXzwwQf069ePqKgoAHr06IHD4eD0009nxIgRle4ndWzRIrj9dkbu2YPH2Hg76XL+esnvyc9oRtbCk8hdHw+WIbLzbpr02Up4/CEAEmJdR4Sk5PbNA255UngSEWk8zOExRaEkOTnZKj8YfP369XTt2jXgx6jtswvrs6q+lvVVhe+B5h743e/gxx8B+LZjb+6+7M/s2ZdI1rcnk7+9FcbpIbrHTpr02Yaj6W8togmxLv435qJgPR0REQlRxpiVlmUll9/eIFuyQF0ujV35M0zTDuSQf8fdWKsXYSwLTj0Vzzvv8cF/4tn4gp2CPU2xR+UTe/4GmvXaiT2ySGsAiohItTTYkCWNW+kzTPtvWs6znzxPk8I83GERMPV13igYynPXwtatEN++CNfg9XgSt5PQMpzRA04teQy1hIqIyPFSyJIGKS3TTdO8LGZ+8DdO37MZC3iry5U82PIZHGO7sm8f9O3rXxHnd79zYrN1Bcp2oSpUiYhIddSrkGVZFqaqsz1KGaE4Bg9qfgzdAz8v4oFPXsHp87KiaU9uaDed1I2nYW1wMnAgPPIInHceVZ48VEREJFD1JmRFRESQkZFBixYtFLSOk2VZZGRkhNzcXDU6Q//u3dC/P6PWrmWzOZmb41/j+73nwVobTU7dzbjH7fxpSOuafgoiIiJHqDchq23btuzatYt9+/YFu5R6LSIigrZt2wa7jDKqOkN/pa1eTz8Njz7KBm9HJjRbwNuHrsTaaxHdfSedLknj8aHtGNRTAUtEROpGvQlZTqeTxMTEYJchtaAqM/RX1Or15KzlXHjTeHZstPOkeZcPzbVE5BsefBAeesgQH98eaF+bT0FEROQI9SZkScNVlRn6y7d69d/4P+5asJg7fOOZy7XERFmM/aNh5Eho1erIY2n+NBERqSu2YBcgMnpAZ1xOe5ltlc1Ldbh1y/i8jH73PQ7M78UZvpUstF/OX/8KO341PPlk5QFr7Nw1pGa6sfht7Nf8Vam18bRERKSRU0uWBF1VFkWOj3VRtHIfbeaGc7/vTaLJIjHpe1pd7Wb8+POPepyqjv0SERGpDoUsCQmBzNC/axe0er8V/151AZHkcUPzmawYGkd4rGHsoO7HPEZVxn6JiIhUl0KWhLyMDJg4weLFyUVYvk7cb14i7KpsPuzakxNjw0u6Fc+ZuPSoLWFVGfslIiJSXQpZErIKC2HKFBg/zkdWFtzKLMY1f5EO374HHTsyqXi/QOfZGj2gc5n9QGsSiohI7dHAdwlJixZBjx4wahSclb2En+jBjIvfpsPeb6FjxzL7Hm2sVWmDeiYwYXB3EmJdGCAh1sWEwd01HktERGqFWrIkpGzc6A9WCxdCp7iDfMItDGQRPPUkjB1b4X2qMtYqkLFfIiIiNUEtWRIS8vPh8cehWzf4+mt4rs+7rElvzQDHEv54+9MkHurBOROXVjjdQmVjqjTWSkREgumYIcsYM90Yk26MWVtq2/XGmHXGGJ8xJvko991ujFljjFltjEmpqaKlYfn6a+jZE/7xD7hpiJdNiQMY9f1N0CSK/r+fxoKWpx51XquqzLMlIiJSVwJpyZoBXFZu21pgMPBlAPfvZ1lWkmVZlYYxaZyysuC+++C888Dthv/MOsCbnyXQ+sdPoVs3Bvz5HbZEl51VVGOtRESkvjjmmCzLsr40xnQot209gDGmdqqSBu/zz2HECEhNhZEj4YlbNhJ9QW/IzYXrr4f33mP72IUV3ldjrUREpD6o7TFZFvCpMWalMebuWj6W1AMeD1x3RzaXXGqxNy+HbvekMKjHR0Sf1d0fsP7yF3j/fTBGY61ERKReq+2QdY5lWb2Ay4H7jDGVrntijLnbGJNijEnZt29fLZclwbBjB3RPLmDO9Biiu+/khOFfc+aeeZx7xzVYRUXw0kvw5JMl+2uslYiI1Ge1GrIsy0or/poOzAPOOMq+0yzLSrYsK7lVRav7Sr02dy4kJcHmDXZaXvUDLS5fw7C1C3hl3lMAPDr0b/4BWqVorJWIiNRntTZPljEmCrBZlpVd/H1/4O+1dTwJTV4vPPIIPPcc9OkDab2+xBHr5oH/zWbU1+9QZLNz85AnSTmxG09VcH+NtRIRkfrqmCHLGDMbuBBoaYzZBfwNOAC8CLQCPjHGrLYsa4AxJh54zbKsgUBrYF7x4HgHMMuyrP/UztOQUJSbCxdd6ea7L1zE9NqO/epttPJ5GLH4bf64/F3yHWFcOfx5trRsR4LGWYmISAMTyNmFN1Vy07wK9k0DBhZ/vxU4vVrVSb0yf1UqkxZvJC3TTQvTlF3v92bvtgiaXbKOJr23szsHRn/5JveueB+3I5wBt09hZ7MTNM5KREQaJC2rIzWi9CLNhXubsPrD3vgKncRd+z2uk/0nMjyy7A3u+W4ObmcEw0a+xi5bLAmxLkYP6KwuQRERaXAUsqRGHF6k2b2lFfs+6oUtoogTbl5OWFw2AH9Z+hp3fT+fPGcEl945leXP3BzkikVERGqXQpbUiLRMN3m/xLFvXm/CWmXT6rrvcUQXAL8FrNwwFxff8TKOdicGuVoREZHap5Alx6X0+Kv4WBdmZxv2zT+dsLgs4oZ8iz3CA8B9y98rCVj97nyFnOZxTND4KxERaQQUsqTKSo+/Ati8MoZ9804nvFUOrUoFrOE/Leahr96iwBHGZbe9SNiJbZmg8VciItJIKGRJlR0efwX4x2DN60VYq2w6Dv+BprFO0jI93Pzrt4xb9CLG6STih5V83a1bkKsWERGpWwpZUmWHF2h2b2lFevEYrLgh35KLh7VjroDPPoMBT4LNBkuXggKWiIg0QrW9dqE0QPGxLgr3NmHfR70Ia5lTMgYrPtYF338Pl10GlgULFsC55wa7XBERkaBQyJIqu71nV/bN6YMtooi4677HHuHB5bTzeFIMnHeefy2dN9+EK64IdqkiIiJBo+5CqZLcXJj6lzY4vD463/k9h1wFxMe6ePiCdlx21ZlQUABPPw233BLsUkVERIJKIUsC5vPBrbfCqlWwYIGNK6/s67/BsuC00yAjA4YPh4cfDm6hIiIiIUAhSwL22GMwdy783//BlVeWuuGqq2D9ejjzTJgxI1jliYiIhBSNyZKAzJwJEybA738PI0eWuuGhh+CTT6BtW/jqq6DVJyIiEmoUsuSY1q71h6uLLoIXXwRjim+YPh2eew6io+HHH8GhhlEREZHDFLLkqNxuuOkmiI2F2bPB6Sy+YdUquOsuf7D65hto3jyodYqIiIQaNT3IUT3yiL8lq+uI1fT9v1TiY12MOTeBqwZe4B8J/8EH/kHvIiIiUoZCllTqk0/83YPNzthOXutUAFIP5tF58ADIzoYHH4TBg4NcpYiISGhSd6FUaM8euO02iDwhm5jz1pdsn7jon3Tat4P1CZ3h+eeDWKGIiEhoU8iSMuavSuXsp5aSeHY6GQe9NL3iB4zDB8A1az9nyJolZEZEM2jIU0GuVEREJLSpu1BKzF+Vyti5a9i7/ETyt8XRvP8awlvmYAEn7d/JpIUv4DU2rrnlWVq2ig12uSIiIiFNLVlSYtLijWTtDefgf7vg6riX6KRfsQCnt4g574zGbvl46IqR7DmhPaMHdA52uSIiIiFNLVlSIvWgmwNLzsDYfbQYsKZkPqzXPxxPs/wcPux2Ed+fcwUTBnRmUM+E4BYrIiIS4hSypET4zvbkb29Fs0vWYo8uAOCm1Ys4b/tqaNeO6376jOtKZiIVERGRo1F3oQCQlQX7lnQlos0hYnruACAhcw9PfPoylt3hXzJHAUtERCRgClkCwOOPQ2aGnScmFdC2uQtjWcybPQa75cP26ivQrl2wSxQREalX1F0o/PADvPQS/OEP8NDNcTzERTBsGGTth8sugzvuCHaJIiIi9Y5asho5r9cfrlq1giefLN740Ufw9tv+9Qj//e+g1iciIlJfqSWrkfvXv+C77/yZKjYWOHQIhgzxj79atsy/ALSIiIhUmVqyGrGMDBg7Fi66CIYOLd542WVQUAB/+xv06BHU+kREROozNVM0Ynf8KZvMQ9FsSPySc5/28qL7B3p98w106uQPWSIiInLcFLIaqdcW7WHBrFZEdduFs2UO2WnZdJs6Fp/Njm3JkmCXJyIiUu+pu7CRevSvPiwg9txNALz1weOEeT28cultmq5BRESkBihkNULr10P6yjbE9NyBo0k+N/+wkNP3/MKmFicyKWlQsMsTERFpEBSyGqHHHgN7mJemZ22hWW4m4z5/FY+xcfOQfxAf6wp2eSIiIg2CxmQ1MPNXpTJp8UbSMt3Ex7oYXW4x5+++g7lz4cZ78ljd1Mtb0/6K0+fliX53kNM8jgkDOgexehERkYZDLVkNyPxVqYydu4bUTDcWkJrpZuzcNcxflQqAZcGYMf6JR6c904RZ/Ei39K1saNme//QfyoTB3csEMhERETl+aslqQCYt3oi7yFtmm7vIy6TFGxnUM4HPPvPPL/rCCxBjz6Pns+PAbqfLj8v5X3x8cIoWERFpoNSS1YCkZbor3e7z+Scebd8efv974IYbID8fHnoIFLBERERqnFqyGpD4WBepFQSt+FgXCxbAypUwYwaEr1wOn3wCrVvDhAl1X6iIiEgjoJasBmT0gM64nPYy21xOOw/178zEiZCYCDcPtWDwYP+Nc+f61ygUERGRGqeWrAbk8KD18mcXNstK4NtvYepUcDw+FvbuhSuugLPPDnLFIiIiDZexLCvYNRwhOTnZSklJCXYZDcbll8MPP8D2FbtxdToRnE7/6tCRkcEuTUREpN4zxqy0LCu5/HZ1FzZwq1fDf/4DDz4IruuvBK/Xf3qhApaIiEitUshq4J55BmJi4N6Ej/zNWV26wN13B7ssERGRBk8hqwHbuhXeew9+f7dF5APD8QEXnP9nzpm4tGSCUhEREakdGvjegD37LDgccNWuxwjLPsTHnc9hR7M2UDwTPKAZ3kVERGqJWrIaqL174Y034NabCjnrg6cptDv48xV/Krn98EzwIiIiUjvUklWPHW0x6H/+EwoKYPSvf8Tp8/LsecMocEaUuX9lM8SLiIhI9Slk1VOHF4M+vFZhaqkuwItOTmDKFBh8SRadlrzKwaimvHTWDUc8Rnysq05rFhERaUyO2V1ojJlujEk3xqwtte16Y8w6Y4zPGHPEvBCl9rvMGLPRGPOLMWZMTRUtR18Mevp0OHQIHvnlLgA2PPkCrrCyedrltDN6QOc6q1dERKSxCWRM1gzgsnLb1gKDgS8ru5Mxxg5MAS4HTgVuMsacenxlSnmVdfWlHnQzZQqc3Xk/fba9D6efzlkPDmfC4O4kxLowQEKsiwmDu2vQu4iISC06ZnehZVlfGmM6lNu2HsAcfd27M4BfLMvaWrzvu8DvgJ+Ps1YppbLFoCP3JbD+Fxgf9Rf/uoRz5gD+swgVqkREROpObZ5dmADsLHV9V/E2qQGVLQYd9UtnWkdlc13uDP9C0CefHJwCRUREGrnaDFkVNXNVulCiMeZuY0yKMSZl3759tVhWwzCoZ8IRXYAj+yax8usI7sp/iTAnMHNmsMsUERFptGrz7MJdwImlrrcF0irb2bKsacA08C8QXYt1NRjluwAffhhsePm9dwr84fcQFRXE6kRERBq32mzJ+h44xRiTaIwJA24EFtTi8Ro1txtef93iGubRNiIDJk8OdkkiIiKNWiBTOMwGVgCdjTG7jDF3GGOuMcbsAs4CPjHGLC7eN94YsxDAsiwPcD+wGFgPvG9Z1rraeiKN3bvvwoEDhvusl/xNWg5NgSYiIhJMxrJCr2cuOTnZSklJCXYZ9YZlQXJPDwU/bmBN9NmYQ5lg04pJIiIidcEYs9KyrCPmDVVzRwPw7bfww48OpjIF8+Q/FLBERERCgEJWA/DSM7k0wcOwFovgjy8HuxwRERGhdge+Sx3Yuxc+mO9kODOJfv4fwS5HREREiqklq56bPjGdQiuOe1vPhVu+CHY5IiIiUkwtWfWYZcHrLxdwAV/QZfrDwS5HRERESlHIqse+fPtXthScyB0tF8DAgcEuR0REREpRyKrHXv/zzzThENdOvTjYpYiIiEg5Cln11KE1v/LhvvMZGvURkddfEexyREREpByFrHpq9g3zcBPJ7WNaB7sUERERqYBCVn20dy/TN5xFd/vPJD86INjViIiISAU0hUM9MX9VKpMWbyQt082492bzPW/z/HUrMCbYlYmIiEhFFLLqgfmrUhk7dw3uIi9N8rL5YXsfwiig2e87BLs0ERERqYS6C+uBSYs34i7yAjD+P9N4m1voGLeGV79dH+TKREREpDIKWfVAWqYbgKiCPPI3x3GAFuy/wFuyXUREREKPQlY9EB/rAmDcZ68wkxHEhu8nosP+ku0iIiISehSy6oHRAzoTZbPouW4bn9Ifq1cGkeF2Rg/oHOzSREREpBIKWfXAoJ4JfHhgKbOtm7GwcfK5GUwY3J1BPROCXZqIiIhUQmcX1hNd3pvBTL6j3/lelj5zdrDLERERkWNQS1Z98O67fJPZmS10ZNgIe7CrERERkQAoZNUHf/kLbzEMV4SPa68NdjEiIiISCIWsULdyJQXbUnnXfjODrrHRpEmwCxIREZFAKGSFugceYCEDOehtyrBhwS5GREREAqWB76Fs71745hveiviE1k3h0kuDXZCIiIgESi1ZoeyBBzhgxfJxUX+GDgWHIrGIiEi9oZAVqoqKYN483g8fRpHXrq5CERGRekZtI6Hqr38Fj4e3Wo3ktOaQlBTsgkRERKQq1JIVguavSiX7+RfZSEeW706k18WHMCbYVYmIiEhVKGSFmPmrUvl0wqvE5OfyVOxowGI5PzB/VWqwSxMREZEqUMgKMZMWb+SBL97CB7zvG0JE+ww8rjwmLd4Y7NJERESkChSyQkzhrlS67NvO/KgryM9qStRpuwBIy3QHuTIRERGpCoWsEPOP/83EAP9o9gjG4SWy0x4A4mNdwS1MREREqkQhK5T4fFyy/muy7FH8uL8vrlP2YAv34nLaGT2gc7CrExERkSpQyAol//oX9oIC5iSNxZcfRvRpaSTEupgwuDuDeiYEuzoRERGpAs2TFUqefhqAT08cRYvtsGtWH5zO4JYkIiIix0ctWaFiwwbYto2cU8/go8UubrgBBSwREZF6TCErVIwaBcD8/lNxu+Hmm4Ncj4iIiFSLQlYo8HhgyRJo2pRZG3vTvj2cdVawixIREZHqUMgKBU89BR4P6dffx6efwtChYNNPRkREpF7Tv/JQMHUq2Gx80PWveL3+kCUiIiL1m0JWsP3wA+zdC3378s6H4XTvDt26BbsoERERqS6FrGB79FEAtt7zDCtWaMC7iIhIQ6GQFUw+H3z+OTRpwuyd5wJw441BrklERERqhEJWME2bBkVFWIOv5Z134LzzoH37YBclIiIiNUEhK5gmTwbgx6FPs369BryLiIg0JApZwZKeDps2QceOzFrSCocDrr8+2EWJiIhITVHICpbiAe++P45k9my47DJo0SLINYmIiEiNUcgKlg8+AKeTL0+9h127dFahiIhIQ6OQFQyffw6HDkG/fjw5JR9bmIdHvl/EOROXMn9VarCrExERkRqgkBUMf/sbAJ8O+RNLF4bhOmUPxukjNdPN2LlrFLREREQagGOGLGPMdGNMujFmbaltzY0xS4wxm4u/Nqvkvl5jzOriy4KaLLzeKiiAb76Bli3549IEfAVOok79LVS5i7xMWrwxiAWKiIhITQikJWsGcFm5bWOAzy3LOgX4vPh6RdyWZSUVX64+/jIbkOeeA68Xhg1j5/ctsUUWENEho8wuaZnuIBUnIiIiNeWYIcuyrC+BA+U2/w6YWfz9TGBQDdfVcL3yChhD5shxuLe0JqprGsZmldklPtYVpOJERESkphzvmKzWlmXtBij+GlfJfhHGmBRjzDfGGAWx7dth507o3p05S5pgeW0077GnzC4up53RAzoHpz4RERGpMbU98L2dZVnJwFDgeWPMyZXtaIy5uziQpezbt6+WywqS8eP9X//0J955Bzp2hMn3tyMh1oUBEmJdTBjcnUE9E4JapoiIiFSf4zjvt9cY08ayrN3GmDZAekU7WZaVVvx1qzHmC6AnsKWSfacB0wCSk5Otivap7wrnzMPY7Jz8XQI7v7C48e5srumVwDW9FKpEREQamuNtyVoADC/+fjjwUfkdjDHNjDHhxd+3BM4Bfj7O49V7y97+hLDsQ6QkdCV7fVuwDN/ZV2u6BhERkQYqkCkcZgMrgM7GmF3GmDuAicClxpjNwKXF1zHGJBtjXiu+a1cgxRjzI7AMmGhZVqMNWc4n/wHAy2deT+7P8YS1OYg3JlvTNYiIiDRQx+wutCzrpkpuuriCfVOAO4u/Xw50r1Z1DYVl0fuXH8hzhPFZzPkUpTel2cXrAE3XICIi0lBpxve68PHHuDyFfJXYk9z1bcH4iOqaBmi6BhERkYZKIasuPP00AK+edzO56+KJ6LAfe1ShpmsQERFpwBSyalthIXz7LcTG0uusG/FmRRJ9WqqmaxAREWngjncKBwnUG2+AxwODB7Plf62IiYEd7/UkMjLYhYmIiEhtUsiqbS+9BEDOyMf44Cy48UYUsERERBoBdRfWpsxMWLcO2rRh7qpEcnNhxIhgFyUiIiJ1QSGrNr3wAlgWDBvGjBlw8slwzjnBLkpERETqgkJWbXrjDQB2DHmYZctg+HAwJsg1iYiISJ1QyKotO3fCjh3QsSNvftICgFtvDXJNIiIiUmcUsmrLc88BYN15FzNnQr9+0L59kGsSERGROqOQVVveew+A//X4A1u2aMC7iIhIY6OQVRt++QX27IFOnZgxJ4boaLj22mAXJSIiInVJIas2vPACAM+3u5jpbxXh6rybJZtSg1yUiIiI1CVNRloL8t+ZTTjwghmKVejE1mk7Y+ceAtAyOiIiIo2EWrJq2tq1RBzMYEdsG/Zs6oyjaR7hJx7AXeRl0uKNwa5ORERE6ohCVk2bOhWAaafcQP6OlkR121UyN1ZapjuIhYmIiEhdUsiqSZZVclbhm57bwEB0j50lN8fHuoJVmYiIiNQxhayalJICBw5woM3J7NnYBVfHvTia5APgctoZPaBzkAsUERGRuqKQVZOmTQPgP2dOxJsXTuK5ezBAQqyLCYO7a9C7iIhII6KzC2uKzwcffADA1B1X0LEjrHktCZstKciFiYiISDCoJaumfPUVHDrEj3GX8r8fXPzhD2DTqysiItJoKQbUlDfeAODl+L8TEaFldERERBo7hayaUFQEc+aQRQxvb+zDTTdB8+bBLkpERESCSSGrJixZAjk5vBVzH7luO/feG+yCREREJNgUsmrCjBlYwFTng/TpA8nJwS5IREREgk0hq7pyc2HBAr7kfH4+cIJasURERATQFA7V99FHUFDAVMeDNIuxGDLEBLsiERERCQFqyaqut98mlXjmeK7G23Ebl7ywlPmrUoNdlYiIiASZQlZ17NuHtXgxz/FnfMYQ3Ws7qZluxs5do6AlIiLSyClkVcf775Pua8kr3EPTLr/ijHUD4C7yMmnxxiAXJyIiIsGkkFUd77zD0+ZhCggn8tztZW5Ky3QHpyYREREJCQpZx2vrVvat2Myr1j10b/0dzuZ5ZW6Oj3UFqTAREREJBQpZx2vWLP6PUbhxEdNvV5mbXE47owd0DlJhIiIiEgoUso6HZZHx5ie8xAMMiV7In0efRUKsCwMkxLqYMLg7g3omBLtKERERCSLNk3U8Vq3i+c0DySGax65Zx2m9rmRQr7bBrkpERERCiFqyjsPB1+fyT/7IdXzAabf2DnY5IiIiEoLUklWJ+atSmbR4I2mZbuJjXYwe0NnfBej18s83Y8miKY+7noPzvwx2qSIiIhKCFLIqMH9VKmPnrsFd5AUomWA0ZccBMuZ8zrs5d3ClbQEtzm4BYWFBrlZERERCkboLKzBp8caSgHWYu8jLO9/8yt4FJ5BJM/7u+xsvRXbRzO4i8v/t3X+sV3Udx/Hnm9slbojA7QIigphDFrUEuwqtYpgIWm1QLas1x9KiOdpqU9NsrWZjWM5yNtdmSVornE0C3drwqpna0gloghoDU+NXIP7GoYC8++OeW9fr98IFbxV5PQAABy9JREFUPPfcr9/nY7s753w43/t9s88++772/XzO50pSTYasGnrbSLRt/RuseOU8Pj9sGafyDzpOnObO7pIkqSZDVg21NhI9sG8QWztOZzybuaTlKh4bO4ldQ0e6s7skSarJkFXDpXMn09Lc9Ja2ffeOZ/ve8Zx/0vVM37mOu08+HXBnd0mSVJsL32vo2ki06+nCY3ePZv3aKVzAjUx93xoGkdw1abo7u0uSpF4Zsnoxf9o45k8bx9690H7aAY6LHfzglN8w8fG/cWP7PF455UMs6drWQZIkqQdD1iEsWQLrHh/EHXyDiU8/DDNncuFdf+TC5uaqS5MkSQOYIesg1q6FxYuTrw5dyWf33wWt74dbbwUDliRJOgQXvvfioYdg9mwYPWwP1772ddi/H5YvhzFjqi5NkiTVAUNWDffcA2edBa2tcP/g2bTxPFx3HcyYUXVpkiSpTjTmdOEll8B9tf/m4MoXZ/KlpxYzachm7tx/PmP/8yBMnw4XXdTPRUqSpHrWp2+yImJpROyMiPXd2lojoiMiNhbHkb28dkFxz8aIWPBOFX5Ujj0W2tre9vO7N87jC0/9hFOHbuKvLecy9tkHYfhwuOMOiKi6akmSVEciMw99U8RMYDfw28z8cNH2U+CFzLwqIi4HRmbmZT1e1wqsBtqBBNYAH83MFw/2fu3t7bl69eoj+f8cttdf75wevO02WLoUPjVmPSt2fIxhrYPhiitg0SIYMqRfapEkSfUnItZkZnvP9j5NF2bmfRExsUfzPGBWcX4zcC9wWY975gIdmflCUUQHcA6wrI91v+NWPLKVlxd+l+f+2cY9++Zw/95Z7M6hDI3dXMAtXPP8xWy74GtMvuZKGDGiqjIlSVKdO5o1WWMycztAZm6PiNE17hkHbO52vaVoq8SKR7byveXr2Lfp+zy1ewqtTbv45PAOZhzzAB9pWcOu4cOZc8YveLV1NEuefo350wxZkiTpyJS98L3WQqaa85MRsRBYCDBhwoRSirl61Qb27HuT1z/zEscNeoDBY1/miWjmCc4Ezvz/jfve5OpVG9zNXZIkHbGj2cJhR0SMBSiOO2vcswUY3+36BGBbrV+WmTdkZntmto8aNeooyurdtpf2ADDkhBd57/EvH3Qte9e9kiRJR+JoQtbtQNfTgguAlTXuWQXMiYiRxdOHc4q2Shw/oqWUeyVJknrq6xYOy4C/A5MjYktEXAhcBZwdERuBs4trIqI9In4NUCx4/zHwcPFzZdci+CpcOncyLc1Nb2lrHhQ0N731K62W5iYunTu5P0uTJEnvMn3awqG/lbmFw4pHtnL1qg1se2kPx49o+V+Y6tnmeixJktQXvW3h0HAhS5Ik6Z3UW8jybxdKkiSVwJAlSZJUAkOWJElSCQxZkiRJJTBkSZIklcCQJUmSVAJDliRJUgkMWZIkSSUwZEmSJJVgQO74HhHPAc+W/DZtwK6S30PlsO/ql31Xn+y3+mXf9Y8TM3NUz8YBGbL6Q0SsrrUFvgY++65+2Xf1yX6rX/ZdtZwulCRJKoEhS5IkqQSNHLJuqLoAHTH7rn7Zd/XJfqtf9l2FGnZNliRJUpka+ZssSZKk0jRkyIqIcyJiQ0RsiojLq65HfRMRz0TEuoh4NCJWV12PehcRSyNiZ0Ss79bWGhEdEbGxOI6sskbV1kvf/SgithZj79GI+HSVNaq2iBgfEX+JiCcj4vGI+HbR7tirSMOFrIhoAq4HzgWmAF+JiCnVVqXDcGZmTvWR5AHvJuCcHm2XA3dn5iTg7uJaA89NvL3vAH5ejL2pmfnnfq5JfbMfuDgzPwjMABYVn2+OvYo0XMgCzgA2Zea/MnMvcAswr+KapHeVzLwPeKFH8zzg5uL8ZmB+vxalPuml71QHMnN7Zq4tzl8FngTG4dirTCOGrHHA5m7XW4o2DXwJ3BkRayJiYdXF6LCNyczt0PlhAIyuuB4dnm9FxGPFdKLTTQNcREwEpgEP4dirTCOGrKjR5iOW9eHjmXkanVO9iyJiZtUFSQ3il8DJwFRgO3BNteXoYCLiGOA24DuZ+UrV9TSyRgxZW4Dx3a5PALZVVIsOQ2ZuK447gT/ROfWr+rEjIsYCFMedFdejPsrMHZn5ZmYeAH6FY2/AiohmOgPW7zNzedHs2KtII4ash4FJEXFSRAwGvgzcXnFNOoSIGBoRw7rOgTnA+oO/SgPM7cCC4nwBsLLCWnQYuj6gC5/DsTcgRUQANwJPZubPuv2TY68iDbkZafH48bVAE7A0MxdXXJIOISI+QOe3VwDvAf5gvw1cEbEMmAW0ATuAHwIrgFuBCcC/gS9mpgusB5he+m4WnVOFCTwDfLNrjY8Gjoj4BHA/sA44UDRfQee6LMdeBRoyZEmSJJWtEacLJUmSSmfIkiRJKoEhS5IkqQSGLEmSpBIYsiRJkkpgyJIkSSqBIUuSJKkEhixJkqQS/BcE0lhlHk9s4QAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 720x360 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "plt.figure(figsize=(10, 5))\n",
    "plt.scatter(ts, ym, label='Data')\n",
    "plt.plot(ts, fopdt(ts, K, tau, theta, y0), color='red', label='FOPDT fit')\n",
    "plt.plot(ts, sopdt(ts, K_2, tau_2, zeta_2, theta_2, y0_2), color='red', label='SOPDT fit')\n",
    "plt.plot(ts, ys + 10, color='blue', label='Original')\n",
    "plt.legend(loc='best')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.7.3"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}